R Notebook

library(tidyverse)

abund <- read.table("../../data/E-WADES/depths_M3.tsv", header=T,
                    row.names=1, sep="\t", check.names=F)
meta_sample <- read.table("../../data/E-WADES/meta_sample.tsv", header=T,
                    row.names=1, sep="\t", check.names=F) %>%
  select(Site, Replicate_Bigger, Date) %>%
  rownames_to_column("sample")
taxon <-  read.table("../../data/E-WADES/taxonomy.tsv", header=T,
                    row.names=1, sep="\t", check.names=F)

all(rownames(abund)==meta_sample$sample) &
all(colnames(abund)==rownames(taxon))

[1] TRUE

Remove replicates from abundance matrix

abund <- abund[meta_sample$Replicate_Bigger,]
meta_sample <- meta_sample[meta_sample$Replicate_Bigger,]

Aggregate to Genus Rank

genus_abund <- abund %>% rownames_to_column("sample") %>%
  pivot_longer(-sample, names_to="Species", values_to="count") %>%
  left_join(taxon%>%rownames_to_column("SpeciesID") %>% select(SpeciesID,Genus) %>%
              rename(Species=SpeciesID), by="Species") %>%
  reframe(count=sum(count), .by=c("sample","Genus")) %>%
  pivot_wider(names_from=Genus, values_from=count) %>%
  column_to_rownames("sample")

genus_rel_abund <- data.frame(genus_abund/rowSums(genus_abund))

subset_genus <- colSums(genus_abund) %>% sort %>% tail(10) %>% names

color_genus <- rownames(qualpalr::qualpal(n = length(subset_genus), 
                                          colorspace = list(h = c(0,360),
                                                            s = c(.25,.6),
                                                            l = c(.45,.85)))$RGB) %>%
  stats::setNames(subset_genus)
color_genus <- c(color_genus, "Other"="#A6A6A6")
color_genus

 Trichococcus Psychrobacter  Giesbergeria     Neisseria Aliarcobacter     Comamonas       Ottowia Acinetobacter 
    "#7CB932"     "#4231B6"     "#B43632"     "#BFD1EF"     "#4879A2"     "#D7A2B0"     "#B939A8"     "#CA8F42" 
Pseudomonas_E    Acidovorax         Other 
    "#42B7A1"     "#EDE6C3"     "#A6A6A6" 

data_composition <- genus_rel_abund %>%
  rownames_to_column("sample") %>%
  pivot_longer(-sample, names_to="genus", values_to="relative") %>%
  mutate(genus = if_else(genus %in% subset_genus, genus, "Other")) %>%
  left_join(meta_sample, by="sample") %>%
  filter(Replicate_Bigger) %>%
  reframe(relative=sum(relative), .by=c("sample","genus","Site")) %>%
  group_by(sample) %>%
  mutate(Pseudomonas_E_rel = sum(relative[genus == "Pseudomonas_E"])) %>%
  ungroup() %>%
  arrange(Pseudomonas_E_rel) %>%
  mutate(sample_fct = factor(sample, levels = unique(sample)))

p.composition <- data_composition %>%
  ggplot(aes(x=sample_fct,y=relative,fill=genus)) +
    geom_bar(stat="identity", color=rgb(0,0,0,.25), linewidth=.25) +
    scale_fill_manual(values=color_genus,
                      labels = function(x) paste0("<i>", x, "</i>")) +
    theme_minimal() +
    theme(legend.position="bottom",
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          legend.text = ggtext::element_markdown()) +
  ylab("Genus\nComposition") +
  labs(fill="Genus") +
  xlab("Samples") +
  ggnewscale::new_scale_fill() +
  geom_tile(data = data_composition,
            aes(x = sample, y = -0.03, fill = Site), 
            height = 0.03, width = 1, inherit.aes = FALSE) 
p.composition

Dimensionality Reduction

# Calculate the distance matrix with bray-curtis and aitchison
dist.bray <- vegan::vegdist(abund, "bray") %>% as.dist
dist.aitc <- abund %>% 
  zCompositions::cmultRepl() %>%
  vegan::vegdist(method = "aitchison") 

No. adjusted imputations:  148 

PCOA

pcoa.bray <- stats::cmdscale(dist.bray, eig=T, add=T)
pcoa.aitc <- stats::cmdscale(dist.aitc, eig=T, add=T)

Plot with Aitchison Distance

centroids.pcoa.aitc <- pcoa.aitc$points %>% as_tibble %>%
  cbind("Site"=meta_sample$Site) %>%
  reframe(V1=mean(V1), 
          V2=mean(V2),
          .by=Site)

varExplained_aitc <- 100 * (pcoa.aitc$eig / sum(pcoa.aitc$eig)) %>%
  round(2)

p.pcoa.aitc <-  pcoa.aitc$points %>% as_tibble %>% 
  cbind("Site"=meta_sample$Site) %>%
  ggplot(aes(x=V1,y=V2, color=Site)) +
  geom_point(size=3) +
  stat_ellipse() +
  geom_point(data=centroids.pcoa.aitc, shape=21, size=5, color="black",
             aes(fill=Site)) +
  theme_minimal() +
  xlab(paste("PCoA1 (~",varExplained_aitc[1],"%)", sep="")) +
  ylab(paste("PCoA2 (~",varExplained_aitc[2],"%)", sep="")) +
  theme(legend.position="bottom") +
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=12))

p.pcoa.aitc

varExplained_bray <- 100 * (pcoa.bray$eig / sum(pcoa.bray$eig)) %>%
  round(2)

p.pcoa.bray <- pcoa.bray$points %>% as_tibble %>% 
  as_tibble() %>% mutate(Pseudomonas_E=genus_rel_abund[, "Pseudomonas_E"]) %>%
  ggplot(aes(x=V1, y=V2, fill=Pseudomonas_E)) +
    geom_point(size=4, shape = 21, color = "black") +
    theme_minimal() +
    xlab(paste("PCoA1 (~",varExplained_bray[1],"%)", sep="")) +
    ylab(paste("PCoA2 (~",varExplained_bray[2],"%)", sep="")) +
    #scale_color_viridis_c(option = "C") +
     scale_fill_gradient(
      low = "#2C3E50",  
      high = color_genus["Pseudomonas_E"]  # Same color associated to Pseudomonas_E genus
    ) +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12)) +
    theme(legend.position="bottom") +
    labs(fill = expression("Relative Composition of "*italic("Pseudomonas_E")))
p.pcoa.bray

Merge the figures

p1 <- p.composition + guides(fill=guide_legend(ncol=6))

p2 <- ggpubr::ggarrange(p.pcoa.bray + theme(legend.margin = margin(t = 20, b = 10)),
                        p.pcoa.aitc, ncol=2, labels=c("B","C"))

pall <- ggpubr::ggarrange(p1, p2, 
                          ncol=1, labels=c("A",""),
                          heights=c(0.4,0.6))
pall

png(filename="E-WADES_pcoa_metrics.png", width=8000, height=6000, res=600)
pall
dev.off()

null device 
          1 

pdf(file="E-WADES_pcoa_metrics.pdf", width=13.3, height=10)

pall
dev.off()

null device 
          1 

Adonis 2

meta_sample <- genus_rel_abund %>% as_tibble(rownames = "sample") %>%
  select(sample, Pseudomonas_E) %>%
  right_join(meta_sample, by = "sample") %>%
  mutate(City = case_when(
    str_detect(Site, "Copenaghen") ~ "Copenaghen",
    TRUE ~ Site
  ))

vegan::adonis2(dist.aitc ~ City + Pseudomonas_E, data = meta_sample, 
               permutations = 10^3, by = "margin")

Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 1000

vegan::adonis2(formula = dist.aitc ~ City + Pseudomonas_E, data = meta_sample, permutations = 10^3, by = "margin")
               Df SumOfSqs      R2      F   Pr(>F)    
City            4   265729 0.32571 32.543 0.000999 ***
Pseudomonas_E   1    35103 0.04303 17.196 0.000999 ***
Residual      229   467480 0.57300                    
Total         234   815846 1.00000                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

vegan::adonis2(dist.bray ~ City + Pseudomonas_E, data = meta_sample, 
               permutations = 10^3, by = "margin")

Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 1000

vegan::adonis2(formula = dist.bray ~ City + Pseudomonas_E, data = meta_sample, permutations = 10^3, by = "margin")
               Df SumOfSqs      R2      F   Pr(>F)    
City            4    7.833 0.16663 15.177 0.000999 ***
Pseudomonas_E   1    6.299 0.13401 48.822 0.000999 ***
Residual      229   29.548 0.62856                    
Total         234   47.009 1.00000                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

UMAP

library(uwot)
set.seed(42)
umap.bray <- uwot::umap(as.dist(dist.bray), min_dist = .5)

set.seed(42)
umap.aitc <- uwot::umap(as.dist(dist.aitc), min_dist = .5)

p.umap.aitc <- umap.aitc %>% as_tibble %>% 
  cbind("Site"=meta_sample$Site) %>%
  ggplot(aes(x=V1,y=V2, fill=Site)) +
  geom_point(size=4, shape = 21, alpha = .75) +
  theme_minimal() +
  xlab("UMAP-1") + ylab("UMAP-2") +
  theme(legend.position="bottom") +
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=12))

p.umap.aitc

p.umap.bray <- umap.bray %>% as_tibble %>% 
    mutate(Pseudomonas_E=genus_rel_abund[, "Pseudomonas_E"]) %>%
    ggplot(aes(x=V1, y=V2, fill=Pseudomonas_E)) +
    geom_point(size=4, shape = 21, color = "black") +
    theme_minimal() +
    xlab("UMAP-1") + ylab("UMAP-2") +
    #scale_color_viridis_c(option = "C") +
     scale_fill_gradient(
      low = "#2C3E50",  
      high = color_genus["Pseudomonas_E"]  # Same color associated to Pseudomonas_E genus
    ) +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12)) +
    theme(legend.position="bottom") +
    labs(fill = expression("Relative Composition of "*italic("Pseudomonas_E")))

p.umap <- ggpubr::ggarrange(
  p.umap.bray + theme(plot.margin = margin(b = 15, t = 5)) + ggtitle("E-WADES Bray-Curtis"), 
  p.umap.aitc + ggtitle("E-WADES Aitchison"), 
  ncol = 2
)
p.umap

png(filename = "../Supplementary/EWADES_umap_metrics.png", width = 3000, height = 1800, res = 300)
p.umap
dev.off()

null device 
          1 

Review

pcoa.JSD <- philentropy::distance(abund, method = "jensen-shannon") %>%
  stats::cmdscale(eig=T, add=T)

Metric: 'jensen-shannon' using unit: 'log'; comparing: 235 vectors.

pcoa.CAN <- vegan::vegdist(abund, method = "canberra") %>%
  stats::cmdscale(eig=T, add=T)
pcoa.HEL <- vegan::vegdist(abund, method = "hellinger") %>%
  stats::cmdscale(eig=T, add=T)

plot.PCOA.pseo <- function(pcoa, meta_sample, genus_abun){
  
  varExplained <- 100 * (pcoa$eig / sum(pcoa$eig)) %>%
    round(2)

  p <- pcoa$points %>% as_tibble %>% 
    as_tibble() %>% mutate(Pseudomonas_E=genus_abun[, "Pseudomonas_E"]) %>%
    ggplot(aes(x=V1, y=V2, color=Pseudomonas_E)) +
    geom_point(size=2) +
    theme_minimal() +
    xlab(paste("PCoA1 (~",varExplained[1],"%)", sep="")) +
    ylab(paste("PCoA2 (~",varExplained[2],"%)", sep="")) +
    scale_fill_gradient(
        low = "#2C3E50",  
        high = color_genus["Pseudomonas_E"]  # Same color associated to Pseudomonas_E genus
      ) +
    theme(legend.position="none") 

  return(p)
}

plot.PCOA.city <- function(pcoa, meta_sample){
  
  varExplained_aitc <- 100 * (pcoa$eig / sum(pcoa$eig)) %>%
    round(2)

  p <-  pcoa$points %>% as_tibble %>% 
    cbind("City" = meta_sample$City) %>%
    ggplot(aes(x=V1,y=V2, color=City)) +
    geom_point(size=2) +
    theme_minimal() +
    xlab(paste("PCoA1 (~",varExplained_aitc[1],"%)", sep="")) +
    ylab(paste("PCoA2 (~",varExplained_aitc[2],"%)", sep="")) +
    theme(legend.position="none") 
}

p.JSD.pseo <- plot.PCOA.pseo(pcoa.JSD, meta_sample, genus_rel_abund)
p.CAN.pseo <- plot.PCOA.pseo(pcoa.CAN, meta_sample, genus_rel_abund)
p.HEL.pseo <- plot.PCOA.pseo(pcoa.HEL, meta_sample, genus_rel_abund)
p.AIT.pseo <- plot.PCOA.pseo(pcoa.aitc, meta_sample, genus_rel_abund)
p.BRA.pseo <- plot.PCOA.pseo(pcoa.bray, meta_sample, genus_rel_abund)

p.JSD.city <- plot.PCOA.city(pcoa.JSD, meta_sample)
p.CAN.city <- plot.PCOA.city(pcoa.CAN, meta_sample)
p.HEL.city <- plot.PCOA.city(pcoa.HEL, meta_sample)
p.AIT.city <- plot.PCOA.city(pcoa.aitc, meta_sample)
p.BRA.city <- plot.PCOA.city(pcoa.bray, meta_sample)

p.pseo <- ggpubr::ggarrange(
  p.JSD.pseo + ggtitle("jensen-shannon pseudomonadaceae"), 
  p.CAN.pseo + ggtitle("canberra pseudomonadaceae"), 
  p.HEL.pseo + ggtitle("hellinger pseudomonadaceae"),
  p.BRA.pseo + ggtitle("bray-curtis pseudomonadaceae"),
  p.AIT.pseo + ggtitle("aitchison pseudomonadaceae"), 
  ncol = 5, legend = "none"
)

p.city <- ggpubr::ggarrange(
  p.JSD.city + ggtitle("jensen-shannon city"), 
  p.CAN.city + ggtitle("canberra city"), 
  p.HEL.city + ggtitle("hellinger city"),
  p.BRA.city + ggtitle("bray-curtis city"),
  p.AIT.city + ggtitle("aitchison city"), 
  ncol = 5, common.legend = T, legend = "none"
)

pall <- ggpubr::ggarrange(p.pseo, p.city, nrow = 2)
pall

png("../Supplementary/EWADES_pcoa_more_metrics.png", width = 6000, height = 2400, res = 300)
pall
dev.off()

null device 
          1